OpenSpaces Data Prep

Sea Surface Trend Preparation for OpenSpaces Presentation

Published

July 12, 2023

Code
gmRi::use_gmri_style_rmd()

Contents:

  1. Progression of Anomalies in NE Atlantic (annual, seasonal)
  2. Representation of Warming Rates
  3. Thermal habitat of BSB
  4. Movement of BSB?

Prepare Monthly SST Data

Looking at it from a thermal preference angle, there is the possibility to look at temperature suitability. This could be summer or winter displays of temperature envelopes that fit their thermal preferences.

Reference on Habitat Suitability https://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0147627&type=printable

Code
# Load Raster
month_sst <- raster::stack(str_c(oisst_path, "monthly_averages/oisst_monthly.nc"), varname = "sst")

# Rename
names(month_sst) <- str_replace_all(str_sub(names(month_sst), 2, -1), "[.]", "-")

# Plot
plot(rotate(month_sst$X1981.09.30), col = temp_pal)

1. Pull September Averages

Code
# Pull Just September:
month_sept <- month_sst[[which(str_sub(names(month_sst), 7, 8) == "09")]]

plot(month_sept$X1981.09.30, col = temp_pal, main = "September 1981")

2. Crop to Our Area

Code
# Cropping Extent  
crop_extent <- make_cropbox(xlims = c(-79, -56), 
                            ylims = c(32, 51))
#st_crs(crop_extent) <- crs(month_sst)

# Crop it
gom_monthly <- mask_nc(month_sept, mask_shape = crop_extent, rotate = T)

# Plot a tester:
plot(gom_monthly$X1981.09.30, col = temp_pal, main = "September 1981")

3. Consolidate Into Decades

Code
# Decades
decades <- list(
  "1990s" = 199,
  "2000s" = 200,
  "2010s" = 201)


# Get the mean temp for each decade
temp_decades <- map(decades, function(decade_id){
  year_layers <- which(str_sub(names(gom_monthly), 2, 4) == decade_id)
  decade_dat <- gom_monthly[[year_layers]]
  decade_mean <- mean(decade_dat, na.rm = T)
  return(decade_mean)
  
  
}) %>% stack()

plot(temp_decades$X1990s, main = "1990's September", col = temp_pal)

4. Cover Data from the Great Lakes

Code
# Make a shape for the area
lakes_cover <- make_cropbox(xlims = c(-80, -75),
                            ylims = c(43, 45))


#lakes_ras <- mask(temp_decades[[1]], lakes_cover, inverse = TRUE)

# Mask that region
temp_decades <- mask(temp_decades, lakes_cover, inverse = TRUE)
plot(temp_decades$X1990s, main = "1990's September", col = temp_pal) 

5. Enhance the Resolution

Code
# Resample to Finer Resolution
# Helps it look better if the grid makes squares
s <- raster(nrow = 1400, 
            ncol = 2000, 
            xmn = extent(crop_extent)[1], 
            xmx = extent(crop_extent)[2], 
            ymn = extent(crop_extent)[3], 
            ymx = extent(crop_extent)[4])

# Resample to higher resolution
temp_hires <- raster::resample(temp_decades, s)

# check one
plot(temp_hires$X1990s, main = "1990's September (Resampled)", col = temp_pal)

6. Mask with Land

Code
# Now that the resolution is sharper we'll want to crop it so that it doesn't look
# So blocky next to the coast:

# Read in shapefile for coasts of US and Canada:
library(rgeoboundaries)

# Canada:
canada_sf <- ne_states("canada", returnclass = "sf")

# USA with MD messed up
# us_sf <- ne_states("united states of america", returnclass = "sf")

# USA with better States
us_sf <- gb_adm1(country = "usa", type = "sscgs")

# Filter to our area
new_england <- us_sf %>% filter(shapeISO %in% str_c("US-", c("ME", "NH", "MA", "CT", "MD", "PA", "NY", "VT", "NC", "SC", "VA", "WV", "DE", "SC", "NJ", "RI")))

# Join them
north_america <- bind_rows(list(mutate(canada_sf, country = "Canada"),
                                mutate(us_sf, country = "United States")))

# North America Zoom Extent
na_bbox <- c(xmin = -82, 
             ymin = 29, 
             xmax = -54, 
             ymax = 52)
north_cropped <- st_crop(north_america, na_bbox)


# # Map it
# ggplot(north_america) + 
#   geom_sf(fill = "gray60", color = "white", size = 0.25) +
#   coord_sf(xlim = c(-82, -53), ylim = c(29, 52), expand = F, crs = 4326) 



# Mask that region
temp_coasts <- mask(temp_hires, north_cropped, inverse = TRUE)
plot(temp_coasts$X1990s, 
     main = "1990's September - coastline smoothed", 
     col = temp_pal, 
     xlim = c(-75, -60),
     ylim = c(40, 48)) 

Decadal Temperature Shifts

The decadal average temperatures for the region are as follows:

1990’s

Code
plot(temp_coasts$X1990s, main = "Average SST 1990's", col = temp_pal)

2000’s

Code
plot(temp_coasts$X2000s, main = "Average SST 2000's", col = temp_pal)

2010’s

Code
plot(temp_coasts$X2010s, main = "Average SST 2010's", col = temp_pal)

Thermal Preference Curves

Knowing that the temperatures globally are changing, it is logical to anticipate species to adjust their movements and behaviors to follow their thermal preferences.

From historical datasets and experiments we are able to determine the range of temperatures that species need and/or prefer in order to live.

Knowing these limits then lets us project onto a map where those conditions exist and where conditions are less-tolerable.

Code
# Get temperature range for study area as a vector
temp_range <- seq(min(values(gom_monthly), na.rm = T),
                  #max(values(gom_monthly), na.rm = T),
                  36,
                  by = .1)




# Function to make preference rasters:
make_pref_ras <- function(in_ras, ras_name, p1, p2, alpha, gamma){
  
  # Make a new raster that we can swap values from  
  pref_ras <- in_ras[[ras_name]]

  # Assign values based on preference curve
  values(pref_ras) <- betaFun(
    values(pref_ras),
    p1 = p1,
    p2 = p2,
    alpha = alpha,
    gamma = gamma)
  
  return(pref_ras)

  
}

Black Seas Bass

Black Sea Bass have a thermal preference (from fall survey data) of 18-25C or 64-77F.

Code
# Use the betaFun(), feed it temperature vector min/max and curve shape
bsb_beta <- betaFun(
  x = temp_range,
        p1 = 18, 
        # p2 = 25, # Lab venture
        p2 = 32, # Expanded thermal preference
        alpha = 1, 
        gamma = 1)


# Show the temp  preference
plot(bsb_beta ~ temp_range, type = "l", main = "Black Sea Bass\nFall Temperature Preference", ylab = "Preference", xlab = "Temperature Range")

Code
# Run thermal preference for each year
hires_list <- unstack(temp_coasts) %>% setNames(names(temp_coasts))
bsb_prefs <- imap(hires_list, make_pref_ras, p1 = 18, p2 = 25, alpha = 1, gamma = 1)

Sea Bass 1990’s

Code
plot(bsb_prefs$X1990s, 
     main = "1990's Black Sea Bass Thermal Preference", 
     col = terrain.colors(14, rev = T), 
     xlim = c(-78, -63), 
     ylim = c(35.5,46))

Code
ggplot() +
  geom_stars(data = st_as_stars(bsb_prefs$X1990s)) +
  geom_sf(data = new_england, fill = "gray90", linewidth = .25) +
    geom_sf(data = canada, fill = "gray90", linewidth = .25) +
    geom_sf(data = greenland, fill = "gray90", linewidth = .25) +
  scale_fill_distiller(
    palette = "YlGnBu", 
    na.value = "transparent", 
    limit = c(0, 1), 
    direction = 1,
    oob = scales::oob_censor,
    breaks = c(0.15, 0.85), 
    labels = c("Not Suitable", "Preferred")) +
    map_theme(
      title = element_text(hjust = 0.5, size = 8),
      axis.text = element_text(size = 6),
      legend.position = "bottom") +
    coord_sf(xlim = c(-75.5, -64), 
             ylim = c(38, 45), 
             expand = F) +
    guides("fill" = guide_colorbar(
      title = "Temperature Preference",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3.5, "in"),
      frame.colour = "black",
      ticks.colour = "black")) +
    labs(title = "Black Sea Bass Habitat\n1990's")

Sea Bass 2000’s

Code
plot(bsb_prefs$X2000s, 
     main = "2000's Black Sea Bass Thermal Preference", 
     col = terrain.colors(14, rev = T), 
     xlim = c(-78, -63), 
     ylim = c(35.5,46))

Sea Bass 2010’s

Code
# plot(bsb_prefs$X2010s, 
#      main = "2010's Black Sea Bass Thermal Preference", 
#      col = terrain.colors(14, rev = T), 
#      xlim = c(-78, -63), 
#      ylim = c(35.5,46))


ggplot() +
  geom_stars(data = st_as_stars(bsb_prefs$X2010s)) +
  geom_sf(data = new_england, fill = "gray90", linewidth = .25) +
    geom_sf(data = canada, fill = "gray90", linewidth = .25) +
    geom_sf(data = greenland, fill = "gray90", linewidth = .25) +
  scale_fill_distiller(
    palette = "YlGnBu", 
    na.value = "transparent", 
    limit = c(0, 1), 
    direction = 1,
    oob = scales::oob_censor,
    breaks = c(0.15, 0.85), 
    labels = c("Not Suitable", "Preferred")) +
    map_theme(
      title = element_text(hjust = 0.5, size = 8),
      axis.text = element_text(size = 6),
      legend.position = "bottom") +
    coord_sf(xlim = c(-75.5, -64), 
             ylim = c(38, 45), 
             expand = F) +
    guides("fill" = guide_colorbar(
      title = "Temperature Preference",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3.5, "in"),
      frame.colour = "black",
      ticks.colour = "black")) +
    labs(title = "Black Sea Bass Habitat\n2010's")

Lobster

Lobster have a thermal preference (lab venture) of 11-22C.

Code
# Use the betaFun(), feed it temperature vector min/max and curve shape
lob_beta <- betaFun(
  x = temp_range,
        p1 = 11, 
        p2 = 22, 
        alpha = 1, 
        gamma = 1)


# Show the temp  preference
plot(lob_beta ~ temp_range, type = "l", main = "Lobster\nFall Temperature Preference", ylab = "Preference", xlab = "Temperature Range")

Code
# Run thermal preference for each year
lob_prefs <- imap(hires_list, make_pref_ras, p1 = 11, p2 = 22, alpha = 1, gamma = 1)

Lobster 1990’s

Code
plot(lob_prefs$X1990s, 
     main = "1990's Lobster Thermal Preference", 
     col = terrain.colors(14, rev = T), 
     xlim = c(-75.25, -60), 
     ylim = c(37.5,48))

Lobster 2000’s

Code
plot(lob_prefs$X2000s, 
     main = "2000's Lobster Thermal Preference", 
     col = terrain.colors(14, rev = T), 
     xlim = c(-75.25, -60), 
     ylim = c(37.5,48))

Lobster 2010’s

Code
# plot(lob_prefs$X2010s, 
#      main = "2010's Lobster Thermal Preference", 
#      col = terrain.colors(14, rev = T), 
#      xlim = c(-75.25, -60), 
#      ylim = c(37.5,48))


ggplot() +
  geom_stars(data = st_as_stars(lob_prefs$X2010s)) +
  geom_sf(data = new_england, fill = "gray90", linewidth = .25) +
    geom_sf(data = canada, fill = "gray90", linewidth = .25) +
    geom_sf(data = greenland, fill = "gray90", linewidth = .25) +
  scale_fill_distiller(
    palette = "YlGnBu", 
    na.value = "transparent", 
    limit = c(0, 1), 
    direction = 1,
    oob = scales::oob_censor,
    breaks = c(0.15, 0.85), 
    labels = c("Not Suitable", "Preferred")) +
    map_theme(
      title = element_text(hjust = 0.5, size = 8),
      axis.text = element_text(size = 6),
      legend.position = "bottom") +
    coord_sf(xlim = c(-75.5, -64), 
             ylim = c(38, 45), 
             expand = F) +
    guides("fill" = guide_colorbar(
      title = "Temperature Preference",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3.5, "in"),
      frame.colour = "black",
      ticks.colour = "black")) +
    labs(title = "Lobster Thermal Habitat\n2010's")

Both Species

1990’s Both

Code
# Leaflet Versions
library(leaflet)
library(leaflet.extras)
library(htmltools)


# Lobster Preference
r <- lob_prefs$X1990s
# BSB Preference
r2 <- bsb_prefs$X1990s


r[r<0.1] <- NA
pal <- colorNumeric(
  #rev(c("#0C2C84", "#41B6C4", "#FFFFCC")),
  #palette = "PuRd",
  palette = "Reds",
  values(r),
  na.color = "transparent")


r2[r2<0.1] <- NA
pal2 <- colorNumeric(
  palette = "Greens",
  #c("#0C2C84", "#41B6C4", "#FFFFCC"), 
  values(r2),
  na.color = "transparent")


# Make a map
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles('Esri.WorldImagery') %>% 
  setView(-68.4, 42, zoom = 6) %>% 
  addRasterImage(
    r, 
    colors = pal, 
    opacity = 0.7) %>%
  addRasterImage(
    r2,
    colors = pal2,
    opacity = 0.7) %>% 
  addLabelOnlyMarkers(
    lng = -68.8, 
    lat = 43.1, 
    label = "Lobster Thermal Habitat", 
    labelOptions = labelOptions(
      noHide = T, 
      style = list(
        "color" = "gray50",
        "font-family" = "serif",
        "font-style" = "bold",
        "font-size" = "12px"))) %>% 
  addLabelOnlyMarkers(
    lng = -72.5, 
    lat = 39.6, 
    label = "Black Sea Bass Thermal Habitat", 
    labelOptions = labelOptions(
      noHide = T, 
      style = list(
        "color" = "gray50",
        "font-family" = "serif",
        "font-style" = "bold",
        "font-size" = "12px"))) 

2010’s Both

Code
# Lobster Preference
r <- lob_prefs$X2010s
# BSB Preference
r2 <- bsb_prefs$X2010s


r[r<0.1] <- NA
pal <- colorNumeric(
  #rev(c("#0C2C84", "#41B6C4", "#FFFFCC")),
  #palette = "PuRd",
  palette = "Reds",
  values(r),
  na.color = "transparent")


r2[r2<0.1] <- NA
pal2 <- colorNumeric(
  palette = "Greens",
  #c("#0C2C84", "#41B6C4", "#FFFFCC"), 
  values(r2),
  na.color = "transparent")


# Make a map
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles('Esri.WorldImagery') %>% 
  setView(-68.4, 42, zoom = 6) %>% 
  addRasterImage(
    r, 
    colors = pal, 
    opacity = 0.7) %>%
  addRasterImage(
    r2,
    colors = pal2,
    opacity = 0.7) %>% 
  addLabelOnlyMarkers(
    lng = -68.8, 
    lat = 43.1, 
    label = "Lobster Thermal Habitat", 
    labelOptions = labelOptions(
      noHide = T, 
      style = list(
        "color" = "gray50",
        "font-family" = "serif",
        "font-style" = "bold",
        "font-size" = "12px"))) %>% 
  addLabelOnlyMarkers(
    lng = -72.5, 
    lat = 39.6, 
    label = "Black Sea Bass Thermal Habitat", 
    labelOptions = labelOptions(
      noHide = T, 
      style = list(
        "color" = "gray50",
        "font-family" = "serif",
        "font-style" = "bold",
        "font-size" = "12px"))) 

Global Warming Rates

Code
# Load warming rates
annual_rates <- raster::stack(str_c(oisst_path, "warming_rates/annual_warming_rates1982to2021.nc"),
                              varname  = "annual_warming_rate")


# Plot Rates
plot(rotate(annual_rates), 
     main = "Annual Warming Rate | 1982-2021", 
     col = terrain.colors(14))

Local Rates

Code
# # Crop to local area:
rates_gom <- mask_nc(ras_obj = annual_rates, mask_shape = crop_extent, rotate = T)

# Mask that region
rates_gom <- mask(rates_gom, lakes_cover, inverse = TRUE)

# Resample to higher resolution
rates_hires <- raster::resample(rates_gom, s)

# Plot Rates
plot(rates_hires, 
     main = "Annual Warming Rate | 1982-2021", 
     col = terrain.colors(14))

Global Warming Ranks

Code
# Plot Ranks:
annual_ranks <- raster::stack(str_c(oisst_path, "warming_rates/annual_warming_rates1982to2021.nc"),
                              varname  = "rate_percentile")
plot(annual_ranks, 
     main = "Warming Rate Percentile | 1982-2021", 
     col = terrain.colors(14))

Local Ranks

Code
# # Crop to local area:
ranks_gom <- mask_nc(ras_obj = annual_ranks, mask_shape = crop_extent, rotate = T)

# Mask that region
ranks_gom <- mask(ranks_gom, lakes_cover, inverse = TRUE)

# Resample to higher resolution
ranks_hires <- raster::resample(ranks_gom, s)

# Plot Ranks
plot(ranks_hires, 
     main = "Warming Rate Percentile | 1982-2021", 
     col = terrain.colors(14))

Bounding the gulf Stream?


Saving Geotiffs

To make the data available to openspaces the requested filetype is a GeoTiff

This stackoverflow question asks details around how those can be saved from R: https://gis.stackexchange.com/questions/257204/saving-geotiff-from-r

More information on GDAL GeoTiff options can be found here: https://gdal.org/drivers/raster/gtiff.html

Code
# Folder location
out_location <- here::here("R/markdown_reports/Openspaces_data")

# List of things to save
tiff_list <- list(
  "sst1990" = temp_coasts$X1990s,
  #"sst2000" = temp_coasts$X2000s,
  "sst2010" = temp_coasts$X2010s,
  "bsb1990" = bsb_prefs$X1990s,
  #"bsb2000" = bsb_prefs$X2000s,
  "bsb2010" = bsb_prefs$X2010s,
  "lob1990" = lob_prefs$X1990s,
  #"lob2000" = lob_prefs$X2000s,
  "lob2010" = lob_prefs$X2010s#,
  # "annualrates" = rates_hires$Annual.Sea.Surface.Temperature.Warming.Rate,
  # "annualranks" = ranks_hires$Annual.Sea.Surface.Temperature.Warming.Rate.Rank
)

# Saving Something
test_save <- tiff_list[1]

Idea 1.: Saving Raster Values as Tiffs

As a starting point for testing, each file was exported as a geotiff file using the values of the data. This was not the desired result, but gave us geo-referenced tiff files that can be formatted using python/GDAL

Code
# # Save a geotiff
# iwalk(tiff_list, function(x, y){
#   writeRaster(x, filename = str_c(out_location, "/", y, ".tif"), options=c('TFW=YES'), overwrite = TRUE)
#   message(str_c("Saving geoTiff: ", y))
# })

Step 2. Add Colormaps to Geotiffs

For the openspaces images we want the rasters to be multi-band color tiffs, not single-band rasters with the data values. The attempt here is to create an encoding for how colors should display

Code
#```{python}

# # Load rasterio
# import rasterio as rio
# import os
# 
# 
# # Check directory:
# os. getcwd()
# 
# # Load one of the fiff files
# input_Dir = 'R/markdown_reports/generated_46.tif'
# src = rasterio.open(input_Dir)
# show(src,cmap="magma")
# 
# 
# # Set Color Levels
# levels = [0, 0.01, 0.02, 0.04, 0.8, 1, 2, 4, 7, 17, 70, 500]
#     clrs = ['#E3E8E81A','#3DFDFF8C', '#3CFAFF', '#6A95FF', '#E61AFE8C', '#FF04FF', '#f4ded9', '#f4dbaa', '#eda14f', '#e87511', '#b55400', '#633b11']
#     
# cmap, norm = colors.from_levels_and_colors(levels, clrs,extend='max')
# 
# show(src,cmap=cmap,norm=norm,interpolation='bilinear')